MODULE SURFWS_INIT_SL_MOD
CONTAINS

SUBROUTINE SURFWS_INIT_SL(KIDIA, KFDIA, KLON, KLEVSN, NCL, PMU0,PSDOR,              &  ! Input
                     &  PTSOIL, PTSKIN,LDLAND, &
                     &  ZDSNTOT, ZSNDEPTH,                  &       
                     &  ZSNPERT,                                           &  ! Input
                     &  ZDSNR,PTSN, PRSN, PSSN, PWSN,PASN,                 &  ! Input
                     &  PTSNWS,PSSNWS,PRSNWS,PWSNWS,                       & ! Output
                     &  PTSNTOP, PTSNBOTTOM, PTSNMIDDLE,                   & ! Output
                     &  PRSNMAX, ZSADEPTH, KLEVSNA, KLEVMID, ZACTDEPTH,    & ! Output
                     &  PTMINCL,PRMINCL,                                   & ! Output
                     &  PTCONSTAVG, PTCONSTSTD,                            & ! Output
                     &  PRCONSTAVG, PRCONSTSTD, PRSNTOP,                   & ! Output
                     &  YDCST, YDSOIL )

USE PARKIND1 , ONLY : JPIM, JPRB
USE YOMHOOK  , ONLY : LHOOK, DR_HOOK, JPHOOK
USE YOS_CST  , ONLY : TCST
USE YOS_SOIL , ONLY : TSOIL

USE ABORT_SURF_MOD

! (C) Copyright 2017- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!**** *SURFWS_INIT_SL* - Snow warm start multi-layer 
!     PURPOSE.
!     --------
!          THIS ROUTINE SETUP PARAMETERS USED IN THE
!          OTHER WARM START ROUTINE BASED ON NSNMLWS VALUE

!**   INTERFACE.
!     ----------
!          *SURFWS_INIT_SL* IS CALLED FROM *SURFWS_CTL*.

!     PARAMETER   DESCRIPTION                                    UNITS
!     ---------   -----------                                    -----
                     
!     INPUT PARAMETERS (INTEGER):
!    *KIDIA*      START POINT
!    *KFDIA*      END POINT
!    *KLON*       Length of arrays
!    *KLEVSN*     Snow vertical levels
!    *KLEVMID*    Snow middle levels (surfws_init)
!    *NCL*        Number of clusters


!     INPUT PARAMETERS (REAL):
!    *ZSNPERT*    snow depth threshold for glaciers
!    *PMU0*       Cosine of solar zenith angle
!    *PSDOR*      sub grid scale orography   (m)
!    *ZDSNTOT*    Total snow depth                            (m)
!    *ZSNDEPTH*   Snow depth of each layer wrt 0              (m)
!    *ZDSNR*      Snow depth per layer (full) wrt 0           (m)

!     INPUT PARAMETERS (LOGICAL):
!    *LDLAND*     LAND/SEA MASK (TRUE/FALSE)

!     INPUT PARAMETERS AT T-1 OR CONSTANT IN TIME (REAL):
!    *PTSOIL*     soil temperature top layer t-1       (K)
!    *PTSKIN*     skin temperature t-1                 (K)
!    *PTSN*       SNOW TEMPERATURE single layer               (K)
!    *PSSN*       SNOW MASS        single layer               (kg m-2)
!    *PRSN*       SNOW DENSITY     single layer               (kg m-3)
!    *PASN*       Snow albedo                                 (K)

!     OUTPUT PARAMETERS (REAL)
!    *PTSNTOP*    Snow temperature top layer                  (K)
!    *PTSNBOTTOM* Snow temperature bottom layer               (K)
!    *PTSNMIDDLE* Snow temperature KLEVMID layer              (K)
!    *PRSNTOP*    Snow density top layer                      (kg m-3)
!    *PRSNMAX*    Snow density MAX allowed                    (kg m-3)
!    *PSADEPTH*   Soil depth top layer                        (m)
!    *PACTDEPTH*  Active snow depth                           (m)
!    *PTCONSTAVG* Constants for temperature exp function
!    *PTCONSTSTD* Constants for temperature energy adj
!    *PRCONSTAVG* Constants for density exp function
!    *PRCONSTSTD* Constants for density energy adj
!    *PRSNTOP*    Snow density top layer

!    OUTPUT PARAMETERS (INTEGER)
!    *KLEVMID*    Snow middle levels 
!    *PTMINCL*    Cluster index for temperature exp function
!    *PRMINCL*    Cluster index for density exp function

!     OUTPUT PARAMETERS (REAL, WARM START):
!    *PTSNWS*        Snow tempeature warm start (initialised)
!    *PRSNWS*        Snow density    warm start (initialised)
!    *PSSNWS*        Snow mass       warm start (initialised)
!    *PWSNWS*        Snow liq water  warm start (initialised)

!     INPUT/OUTPUT PARAMETERS
!    *KLEVSNA*    Snow vertical Active levels

!     METHOD.
!     -------
!     Values of exp functions are pre-computed with k-cluster 
!     algorithm, see Arduini et al. 2019. 
!     Snow density top layer is computed with linear regression,
!     with pre-computed parameters.

!     EXTERNALS.
!     ----------
!          NONE.

!     REFERENCE.
!     ----------
!          Arduini et al. (2019)

!     Modifications:
!     Original   G. Arduini      ECMWF     28/07/2017

!     ------------------------------------------------------------------

IMPLICIT NONE

! Declaration of arguments 
INTEGER(KIND=JPIM), INTENT(IN) :: KIDIA
INTEGER(KIND=JPIM), INTENT(IN) :: KFDIA
INTEGER(KIND=JPIM), INTENT(IN) :: KLON
INTEGER(KIND=JPIM), INTENT(IN) :: KLEVSN
INTEGER(KIND=JPIM), INTENT(IN) :: NCL
LOGICAL,            INTENT(IN) :: LDLAND(:)

REAL(KIND=JPRB), INTENT(IN)    :: ZDSNTOT(:), PTSN(:), PRSN(:), PSSN(:), PWSN(:)
REAL(KIND=JPRB), INTENT(IN)    :: ZSNDEPTH(:,:)
REAL(KIND=JPRB), INTENT(IN)    :: ZDSNR(:,:)
REAL(KIND=JPRB), INTENT(IN)    :: PSDOR(:)

REAL(KIND=JPRB), INTENT(IN)    :: PMU0(:)
REAL(KIND=JPRB), INTENT(IN)    :: PTSOIL(:)
REAL(KIND=JPRB), INTENT(IN)    :: PTSKIN(:)
REAL(KIND=JPRB), INTENT(IN)    :: PASN(:)

TYPE(TCST)     , INTENT(IN) :: YDCST
TYPE(TSOIL)    , INTENT(IN) :: YDSOIL

! Output Variables
REAL(KIND=JPRB), INTENT(OUT)    :: PTSNWS(:,:)
REAL(KIND=JPRB), INTENT(OUT)    :: PSSNWS(:,:)
REAL(KIND=JPRB), INTENT(OUT)    :: PRSNWS(:,:)
REAL(KIND=JPRB), INTENT(OUT)    :: PWSNWS(:,:)

REAL(KIND=JPRB), INTENT(OUT)    :: PTSNTOP(:), PTSNBOTTOM(:), PTSNMIDDLE(:)
REAL(KIND=JPRB), INTENT(OUT)    :: PRSNTOP(:)

REAL(KIND=JPRB), INTENT(OUT)    :: PRSNMAX(:)
REAL(KIND=JPRB), INTENT(OUT)    :: ZSADEPTH(:)

! Active number of layers:
INTEGER(KIND=JPIM),INTENT(INOUT)  :: KLEVSNA(:)

! Active depth layer (glaciers):
REAL(KIND=JPRB)    :: ZACTDEPTH(:)

! Cluster values:
INTEGER(KIND=JPIM),            INTENT(OUT)  :: KLEVMID(:)
INTEGER(KIND=JPIM),            INTENT(OUT)  :: PTMINCL(:), PRMINCL(:) 
REAL(KIND=JPRB),               INTENT(OUT)  :: PRCONSTAVG(:,:), PRCONSTSTD(:,:)
REAL(KIND=JPRB),               INTENT(OUT)  :: PTCONSTAVG(:,:), PTCONSTSTD(:,:)

REAL(KIND=JPRB)                 :: ZSNPERT

! Local variables:
INTEGER(KIND=JPIM)              :: JL,JK,KNACC
INTEGER(KIND=JPIM)              :: KCL
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTDIST, ZRDIST
REAL(KIND=JPRB)                 :: ZFEATA, ZFEATB, ZFEATC 

! Look-up tables for warm-start
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTCENTRADAY2, ZTCENTRBDAY2 ,ZTCENTRCDAY2 ,ZTCENTRADAY3 ,ZTCENTRBDAY3 ,ZTCENTRCDAY3 ,&
                                    & ZTCENTRADAY4 ,ZTCENTRBDAY4 ,ZTCENTRCDAY4 ,ZTCENTRADAY5 ,ZTCENTRBDAY5 ,ZTCENTRCDAY5 ,&
                                    & ZTCENTRADAY5M ,ZTCENTRBDAY5M ,ZTCENTRCDAY5M,&
                                    & ZTCENTRADAY5G ,ZTCENTRBDAY5G ,ZTCENTRCDAY5G
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTCENTRANIGHT2 ,ZTCENTRBNIGHT2 ,ZTCENTRCNIGHT2 ,ZTCENTRANIGHT3 ,&
                                    & ZTCENTRBNIGHT3 ,ZTCENTRCNIGHT3 ,ZTCENTRANIGHT4 ,ZTCENTRBNIGHT4 ,&
                                    & ZTCENTRCNIGHT4 ,ZTCENTRANIGHT5 ,ZTCENTRBNIGHT5 ,ZTCENTRCNIGHT5 ,&
                                    & ZTCENTRANIGHT5M ,ZTCENTRBNIGHT5M ,ZTCENTRCNIGHT5M, &
                                    & ZTCENTRANIGHT5G ,ZTCENTRBNIGHT5G ,ZTCENTRCNIGHT5G
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTCONSTAVGDAY2 ,ZTCONSTSTDDAY2 ,ZTCONSTAVGDAY3 ,ZTCONSTSTDDAY3 ,&
                                    & ZTCONSTAVGDAY4 ,ZTCONSTSTDDAY4 ,ZTCONSTAVGDAY5 ,ZTCONSTSTDDAY5 ,&
                                    & ZTCONSTAVGDAY5M ,ZTCONSTSTDDAY5M, &
                                    & ZTCONSTAVGDAY5G ,ZTCONSTSTDDAY5G
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTCONSTAVGNIGHT2 ,ZTCONSTSTDNIGHT2 ,ZTCONSTAVGNIGHT3 ,ZTCONSTSTDNIGHT3,&
                                    & ZTCONSTAVGNIGHT4 , ZTCONSTSTDNIGHT4 ,ZTCONSTAVGNIGHT5 ,ZTCONSTSTDNIGHT5,&
                                    & ZTCONSTAVGNIGHT5M ,ZTCONSTSTDNIGHT5M, &
                                    & ZTCONSTAVGNIGHT5G ,ZTCONSTSTDNIGHT5G

REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTCONSTAVG, ZRCONSTAVG, ZTCONSTSTD, ZRCONSTSTD
REAL(KIND=JPRB), DIMENSION(NCL)  ::  ZTCENTRA, ZTCENTRB, ZTCENTRC
REAL(KIND=JPRB), DIMENSION(NCL)  ::  ZRCENTRA, ZRCENTRB, ZRCENTRC


REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTMLRIDAY2 ,ZTMLRADAY2 ,ZTMLRBDAY2 ,ZTMLRCDAY2 ,ZTMLRDDAY2 ,ZTMLREDAY2 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTMLRINIGHT2 ,ZTMLRANIGHT2 ,ZTMLRBNIGHT2 ,ZTMLRCNIGHT2 ,ZTMLRDNIGHT2 ,ZTMLRENIGHT2
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTMLRIDAY3 ,ZTMLRADAY3 ,ZTMLRBDAY3 ,ZTMLRCDAY3 ,ZTMLRDDAY3 ,ZTMLREDAY3
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTMLRINIGHT3 ,ZTMLRANIGHT3 ,ZTMLRBNIGHT3 ,ZTMLRCNIGHT3 ,ZTMLRDNIGHT3 ,ZTMLRENIGHT3
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCENTRA2 ,ZRCENTRB2 ,ZRCENTRC2 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCENTRA3 ,ZRCENTRB3 ,ZRCENTRC3 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCENTRA4 ,ZRCENTRB4 ,ZRCENTRC4
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCENTRA5 ,ZRCENTRB5 ,ZRCENTRC5
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCENTRA5M ,ZRCENTRB5M ,ZRCENTRC5M
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCONSTAVG2 ,ZRCONSTSTD2 ,ZRCONSTAVG3 ,ZRCONSTSTD3 ,&
                                    ZRCONSTAVG4 ,ZRCONSTSTD4 ,ZRCONSTAVG5 ,ZRCONSTSTD5 ,ZRCONSTAVG5M ,ZRCONSTSTD5M 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI2 ,ZRMLRA2 ,ZRMLRB2 ,ZRMLRC2 ,ZRMLRD2 ,ZRMLRE2 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI3 ,ZRMLRA3 ,ZRMLRB3 ,ZRMLRC3 ,ZRMLRD3 ,ZRMLRE3 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI4 ,ZRMLRA4 ,ZRMLRB4 ,ZRMLRC4 ,ZRMLRD4 ,ZRMLRE4 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI5 ,ZRMLRA5 ,ZRMLRB5 ,ZRMLRC5 ,ZRMLRD5 ,ZRMLRE5 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI5M ,ZRMLRA5M ,ZRMLRB5M ,ZRMLRC5M ,ZRMLRD5M ,ZRMLRE5M

REAL(KIND=JPRB), DIMENSION(NCL)  ::  ZTMLRA, ZTMLRB, ZTMLRC, ZTMLRD, ZTMLRE, ZTMLRI
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI, ZRMLRA, ZRMLRB, ZRMLRC, ZRMLRD, ZRMLRE
REAL(KIND=JPRB)                  :: ZP1, ZP2, ZP3, ZP4, ZP5
REAL(KIND=JPRB)                  :: ZSDORTHR
REAL(KIND=JPRB)                  :: ZEPSILON

REAL(KIND=JPHOOK)                  :: ZHOOK_HANDLE


!    -----------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('SURFWS_INIT_SL_MOD:SURFWS_INIT_SL',0,ZHOOK_HANDLE)

!    -----------------------------------------------------------------

ASSOCIATE(RTT=>YDCST%RTT,RLMLT=>YDCST%RLMLT, RPI=>YDCST%RPI,    &
        & RLWCSWEA=>YDSOIL%RLWCSWEA, RLWCSWEB=>YDSOIL%RLWCSWEB, &
        & RLWCSWEC=>YDSOIL%RLWCSWEC, RTEMPAMP=>YDSOIL%RTEMPAMP, &
        & RDSNMAX=>YDSOIL%RDSNMAX, RHOMINSND=>YDSOIL%RHOMINSND, &
        & RHOMAXSN_NEW=>YDSOIL%RHOMAXSN_NEW, RDAT=>YDSOIL%RDAT  )
! mid-point first soil level
ZSADEPTH(KIDIA:KFDIA)=0.5_JPRB*RDAT(1)
ZSDORTHR=50._JPRB

ZEPSILON  = 10E4*EPSILON(ZEPSILON)

! 0.1 Define centroids 
! A: soT-skt B: soT-snT ; C: rsn/rsnmin
    ZTCENTRADAY2(1:NCL)=(/ -7.981,6.879,17.362,1.512,-1.438,26.081,11.509,-15.51,4.07,-3.196,18.445,0.175,8.175,71.994,13.221,22.386,-5.947,-3.131,-11.46,15.301,-0.784,9.829,5.342,19.493,2.636,3.463,-5.436  /)
    ZTCENTRBDAY2(1:NCL)=(/ -2.174,5.977,15.554,1.754,-0.186,23.244,10.122,-2.092,3.92,-1.411,3.092,0.667,7.323,5.033,11.747,19.765,-0.078,-0.083,-2.226,13.6,0.017,8.636,4.963,17.333,3.115,1.448,-1.977  /)
    ZTCENTRCDAY2(1:NCL)=(/ 3.031,2.156,1.977,2.217,2.339,1.847,2.044,3.214,2.096,2.527,1.972,2.254,2.077,1.595,2.038,1.901,4.744,4.868,3.247,2.015,4.844,2.097,2.111,1.978,2.172,4.005,2.655  /)
    ZTCENTRADAY3(1:NCL)=(/ -4.153,13.557,3.599,24.981,0.952,9.803,19.651,-9.364,6.402,-2.604,33.167,28.508,5.363,15.627,8.1,-2.137,2.216,-0.255,-6.922,-5.302,5.008,17.458,-11.989,-15.296,-0.645,22.091,11.631  /)
    ZTCENTRBDAY3(1:NCL)=(/ -1.934,12.535,3.864,23.003,1.413,9.253,18.145,-2.127,6.18,-0.102,30.616,26.221,2.078,14.488,7.523,-0.897,2.728,0.383,-2.067,-0.036,4.89,16.299,-2.351,-2.312,0.401,20.35,10.857  /)
    ZTCENTRCDAY3(1:NCL)=(/ 2.62,2.06,2.189,1.881,2.372,2.087,1.926,3.183,2.152,4.762,1.805,1.818,4.839,1.94,2.189,2.546,2.333,4.843,2.996,4.55,2.189,1.911,3.205,3.136,2.255,1.901,2.092  /)
    ZTCENTRADAY4(1:NCL)=(/ -5.648,9.287,24.102,5.496,16.175,4.544,19.689,-0.478,-4.78,12.576,-10.527,7.71,27.07,1.639,-2.676,-2.336,2.992,-0.323,-14.372,14.384,21.735,18.002,10.921,-0.459,31.438,-7.541,6.187  /)
    ZTCENTRBDAY4(1:NCL)=(/ 0.337,9.207,22.751,1.648,15.55,5.008,18.86,2.65,-1.71,12.178,-2.189,7.727,25.374,2.078,-0.013,-0.909,3.797,0.499,-2.24,13.834,20.637,17.201,10.608,0.207,29.333,-2.155,6.239  /)
    ZTCENTRCDAY4(1:NCL)=(/ 4.804,2.239,1.91,4.57,2.027,2.345,1.92,2.732,2.664,2.149,3.235,2.32,1.881,2.493,4.417,2.451,2.407,2.368,3.217,2.107,1.908,1.998,2.208,4.815,1.874,2.91,2.284  /)
    ZTCENTRADAY5(1:NCL)=(/ 9.309,-2.459,14.928,24.975,4.097,8.582,-8.612,19.429,3.848,6.199,0.459,35.13,12.709,-5.715,-5.041,-1.97,-2.501,5.111,2.558,-12.878,17.164,22.022,7.047,29.34,11.267,0.155,0.644  /)
    ZTCENTRBDAY5(1:NCL)=(/ 8.768,0.513,15.087,24.299,4.718,11.449,-1.685,19.338,7.328,6.139,1.201,34.263,13.263,0.615,-1.151,-0.811,1.641,1.252,2.983,-1.835,17.202,21.698,8.243,28.409,10.954,0.619,4.296  /)
    ZTCENTRCDAY5(1:NCL)=(/ 2.333,5.066,2.134,1.963,2.436,2.35,3.064,2.039,2.585,2.347,2.557,1.87,2.227,5.098,2.881,3.452,2.805,4.479,2.519,3.334,2.102,2.008,2.343,1.965,2.259,4.834,2.692  /)
    ZTCENTRADAY5M(1:NCL)=(/ 12.87,1.494,9.417,5.791,20.58,-1.894,30.502,2.093,9.085,7.134,4.279,9.926,24.749,-4.142,12.111,7.772,13.06,5.954,-0.078,19.413,4.284,15.751,3.009,-0.007,-7.313,16.161,-1.295  /)
    ZTCENTRBDAY5M(1:NCL)=(/ 15.542,1.625,9.898,4.91,20.524,0.846,29.919,8.752,13.191,7.863,1.463,6.204,24.105,0.692,12.272,2.123,9.361,10.676,4.755,16.197,7.016,13.125,4.121,1.034,0.543,17.787,-0.336  /)
    ZTCENTRCDAY5M(1:NCL)=(/ 2.739,4.043,3.025,3.236,2.416,4.849,2.11,3.65,3.022,3.143,4.446,3.576,2.298,4.909,2.82,4.424,3.404,3.272,3.684,2.848,3.376,3.031,3.35,4.651,4.87,2.581,3.552  /)


    ZTCENTRADAY5G(1:NCL)=(/ -4.295,-14.778,6.153,-7.977,10.353,-19.544,-11.999,3.32,-7.696,-10.912,0.809,-4.2,2.475,-15.163,-17.261,5.282,-9.299,8.696,-0.276,-1.255,13.841,-2.726,-22.621,-5.892,-13.383,-6.829,-10.534  /)
    ZTCENTRBDAY5G(1:NCL)=(/ -4.117,-11.803,3.933,-7.408,9.916,-15.274,-12.677,1.222,-4.602,-6.728,-1.467,-7.247,4.589,-14.947,-12.848,7.35,-11.898,7.074,1.797,-4.623,12.157,-0.732,-19.885,-1.56,-9.431,-9.728,-9.469  /)
    ZTCENTRCDAY5G(1:NCL)=(/ 3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.748,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75  /)




    ZTCENTRANIGHT2(1:NCL)=(/ 14.448,3.338,21.303,8.68,-1.001,11.335,7.264,26.333,17.607,0.539,7.359,10.011,-2.748,-5.145,0.502,3.885,12.891,23.431,1.871,49.874,24.46,15.971,4.745,19.351,6.011,-8.798,31.407  /)
    ZTCENTRBNIGHT2(1:NCL)=(/ 13.016,3.226,19.368,7.792,-0.568,10.273,3.342,23.944,16.001,0.715,6.562,9.017,-1.71,-2.102,0.282,1.22,11.667,21.372,1.893,4.817,4.347,14.481,4.389,17.532,5.472,-2.009,29.272  /)
    ZTCENTRCNIGHT2(1:NCL)=(/ 2.03,2.146,1.935,2.125,2.397,2.048,4.132,1.898,1.991,2.165,2.132,2.087,2.593,2.581,4.685,4.766,2.043,1.889,2.153,1.788,1.798,2.019,2.123,1.982,2.112,2.929,1.645  /)
    ZTCENTRANIGHT3(1:NCL)=(/ 1.657,17.38,33.759,8.507,24.649,11.407,-4.348,4.544,-1.798,14.531,5.848,20.67,-0.052,6.675,8.815,0.372,29.796,22.414,10.017,12.912,7.044,16.011,3.264,-8.039,3.148,18.848,27.037  /)
    ZTCENTRBNIGHT3(1:NCL)=(/ 1.688,16.102,31.881,7.943,23.011,10.527,-1.88,4.318,-1.349,13.417,5.366,19.205,0.259,1.434,5.422,0.1,28.013,21.002,9.18,11.939,6.625,14.773,0.84,-1.832,3.058,17.501,25.185  /)
    ZTCENTRCNIGHT3(1:NCL)=(/ 2.247,1.958,1.776,2.138,1.863,2.087,2.78,2.256,2.689,2.0,2.196,1.928,2.256,5.186,3.193,4.622,1.787,1.899,2.13,2.062,2.135,1.972,4.722,2.702,2.254,1.971,1.844  /)
    ZTCENTRANIGHT4(1:NCL)=(/ -2.556,16.089,8.436,22.12,6.305,28.095,11.608,5.586,2.55,35.31,18.984,0.139,-0.443,13.155,6.902,-6.589,23.799,20.517,17.551,4.209,10.112,2.264,14.519,30.57,25.826,9.595,0.83  /)
    ZTCENTRBNIGHT4(1:NCL)=(/ -1.606,15.038,7.983,20.926,1.339,26.63,10.707,5.127,2.58,33.489,17.909,0.259,-0.462,12.136,6.473,-2.189,22.505,19.362,16.437,4.006,9.323,0.199,13.56,28.953,24.434,5.885,1.089  /)
    ZTCENTRCNIGHT4(1:NCL)=(/ 2.495,2.111,2.225,1.934,4.723,1.853,2.223,2.321,2.335,1.788,1.989,4.917,2.363,2.181,2.302,2.602,1.902,1.972,2.087,2.341,2.254,3.81,2.129,1.858,1.864,2.885,2.244  /)
    ZTCENTRANIGHT5(1:NCL)=(/ 0.512,19.694,11.149,27.953,4.797,16.347,23.396,12.846,8.807,1.476,-1.716,36.275,21.422,9.377,31.345,7.277,13.747,7.776,0.813,18.013,3.145,-6.314,6.215,14.742,25.446,11.245,4.014  /)
    ZTCENTRBNIGHT5(1:NCL)=(/ -0.192,18.346,10.503,26.399,4.366,15.144,21.965,12.348,5.269,1.612,-1.147,34.234,20.181,8.897,29.618,1.439,10.649,7.344,0.615,16.761,3.028,-1.72,5.839,13.564,23.993,7.899,0.719  /)
    ZTCENTRCNIGHT5(1:NCL)=(/ 2.811,2.07,2.208,1.958,2.428,2.146,2.028,2.147,2.753,2.46,2.887,1.88,2.012,2.256,1.904,4.599,2.516,2.297,5.111,2.112,2.488,3.075,2.336,2.2,1.98,2.61,4.519  /)
    ZTCENTRANIGHT5M(1:NCL)=(/ 5.949,17.746,12.328,30.563,-0.649,3.874,10.352,23.518,13.354,16.227,19.184,8.179,15.303,8.78,6.215,9.565,15.08,26.622,19.617,5.292,10.839,12.436,22.168,3.357,1.49,7.573,35.391  /)
    ZTCENTRBNIGHT5M(1:NCL)=(/ 3.767,13.235,10.385,28.564,-0.008,0.974,7.96,21.502,8.093,16.087,18.632,5.429,10.655,1.952,1.216,10.513,13.173,24.535,15.672,6.114,5.138,13.034,18.281,3.402,1.02,7.998,33.1  /)
    ZTCENTRCNIGHT5M(1:NCL)=(/ 3.417,3.167,3.051,2.126,4.0,4.595,3.358,2.306,3.661,2.514,2.381,3.45,3.4,4.466,4.654,2.778,2.877,2.262,2.958,3.038,3.783,2.611,2.71,3.182,4.378,2.874,2.03  /)


    ZTCENTRANIGHT5G(1:NCL)=(/ 8.646,-4.116,17.832,0.864,13.553,3.306,-13.676,4.103,-3.558,15.345,7.959,10.417,6.273,-6.906,0.014,0.08,9.283,-3.586,12.804,-8.389,10.967,-0.992,6.082,2.21,20.907,14.257,3.953  /)
    ZTCENTRBNIGHT5G(1:NCL)=(/ 5.237,-6.464,15.365,1.736,9.252,4.252,-10.678,8.443,-2.222,12.488,8.495,7.322,2.912,-3.087,-3.215,5.289,12.155,1.65,11.715,-8.374,9.672,-0.374,5.802,-0.713,18.464,14.556,1.457  /)
    ZTCENTRCNIGHT5G(1:NCL)=(/ 3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.75,3.749,3.75,3.75,3.75,3.75,3.75  /)


  ZTCONSTAVGDAY2(1:NCL)=(/ -0.1,12.338,9.796,3.439,-0.908,8.98,13.406,-0.1,7.781,-1.186,9.49,0.043,14.778,-0.1,13.568,9.597,-0.118,-0.137,-0.1,12.222,-0.123,21.834,11.454,10.18,7.013,16.779,-0.248  /)
    ZTCONSTSTDDAY2(1:NCL)=(/ 0.0,31.485,15.942,19.744,5.903,11.726,30.836,0.0,27.289,7.187,16.031,5.249,34.751,0.0,24.085,16.623,1.053,0.882,0.0,22.027,4.783,153.431,32.051,18.141,26.104,53.017,2.954  /)
    ZTCONSTAVGDAY3(1:NCL)=(/ -4.043,14.62,6.372,8.459,2.88,13.769,9.757,-3.945,9.082,-0.141,8.121,6.746,2.181,11.458,9.732,-1.471,5.26,0.136,-5.133,-0.187,5.948,11.586,-4.478,-6.176,0.912,10.36,13.14  /)
    ZTCONSTSTDDAY3(1:NCL)=(/ 12.619,27.31,19.541,12.408,12.687,26.803,15.608,12.446,22.417,3.83,8.525,8.04,10.86,19.807,22.091,7.927,17.003,5.77,14.509,4.019,17.871,21.25,13.108,14.858,7.854,18.055,25.78  /)
    ZTCONSTAVGDAY4(1:NCL)=(/ -0.086,22.317,8.511,0.491,13.692,15.139,9.807,24.705,2.16,18.568,7.08,19.454,10.735,12.673,0.253,1.235,14.669,2.469,8.685,17.698,10.018,13.53,19.883,-0.059,7.648,5.41,16.776  /)
    ZTCONSTSTDDAY4(1:NCL)=(/ 0.116,47.609,9.19,2.629,27.705,39.647,17.701,51.474,5.732,38.654,9.539,44.649,29.576,36.143,6.883,5.228,40.044,17.5,10.524,39.262,16.847,29.99,42.521,1.322,6.937,8.117,41.572  /)
    ZTCONSTAVGDAY5(1:NCL)=(/ 11.932,-0.09,14.309,11.823,15.219,21.996,1.861,11.457,28.074,12.574,7.57,8.589,14.535,-0.1,0.366,0.188,8.338,0.834,13.295,2.869,14.27,10.054,21.301,11.737,13.563,-0.305,42.687  /)
    ZTCONSTSTDDAY5(1:NCL)=(/ 29.946,0.107,30.986,23.163,38.433,43.616,4.146,20.496,51.594,33.045,26.781,6.672,33.005,0.024,2.738,3.151,29.812,8.613,34.957,4.587,31.829,16.795,43.235,24.499,31.794,2.631,58.68  /)
    ZTCONSTAVGDAY5M(1:NCL)=(/ 4.092,1.574,2.896,3.172,4.502,-0.012,5.956,2.183,2.791,4.609,1.322,2.764,5.717,-0.084,3.426,1.417,2.179,3.338,9.068,2.54,6.004,2.355,4.092,0.366,-0.088,4.776,0.684  /)
    ZTCONSTSTDDAY5M(1:NCL)=(/ 13.812,9.856,9.449,11.799,6.125,2.885,3.391,17.196,9.869,17.006,6.435,6.786,10.411,0.324,11.282,5.916,2.456,15.442,28.798,1.85,22.433,1.808,15.973,5.369,0.499,13.048,7.232  /)

    ZTCONSTAVGDAY5G(1:NCL)=(/ 3.254,1.32,3.923,1.652,2.476,1.519,1.123,1.527,3.24,2.529,1.588,1.244,3.418,1.337,1.47,2.635,1.051,3.436,2.799,2.023,2.339,3.708,1.839,4.5,1.611,1.133,1.235  /)
    ZTCONSTSTDDAY5G(1:NCL)=(/ 5.291,1.885,5.595,2.798,2.462,1.409,0.549,4.186,5.098,4.11,4.087,1.699,6.228,0.639,2.163,2.492,0.404,4.209,7.69,4.029,1.605,6.241,0.851,6.205,2.614,0.963,1.527  /)

   ZTCONSTAVGNIGHT2(1:NCL)=(/ 8.231,4.698,7.679,8.496,-0.056,8.843,75.161,6.501,8.021,0.483,6.886,9.727,-0.078,-0.065,1.327,9.684,8.599,7.576,2.237,71.547,14.194,8.302,4.875,7.956,6.404,-0.1,6.35  /)
    ZTCONSTSTDNIGHT2(1:NCL)=(/ 11.938,18.633,9.712,20.38,5.064,14.412,361.395,9.117,10.29,7.073,16.474,28.854,4.928,1.749,12.604,44.803,12.847,9.767,14.685,21.136,20.12,10.586,15.749,10.109,18.877,0.0,8.21  /)
    ZTCONSTAVGNIGHT3(1:NCL)=(/ 2.084,8.251,7.999,8.004,7.885,8.709,-3.302,5.534,-2.179,9.188,6.186,8.955,0.256,-0.017,3.778,-0.32,8.258,7.992,8.833,9.596,7.594,8.039,1.002,-4.867,4.947,9.572,6.438  /)
    ZTCONSTSTDNIGHT3(1:NCL)=(/ 10.173,13.827,8.578,17.67,10.764,17.282,11.497,15.952,9.37,17.872,16.072,13.863,6.472,0.379,29.121,2.766,9.291,10.053,18.707,19.095,18.137,14.253,7.818,14.089,15.659,13.54,8.682  /)
    ZTCONSTAVGNIGHT4(1:NCL)=(/ 3.468,7.444,6.203,6.6,-0.078,6.596,6.941,4.961,5.394,7.086,6.385,-0.031,-0.669,7.355,6.255,5.464,6.741,6.762,6.651,4.599,6.442,0.109,7.052,6.72,6.44,0.195,4.174  /)
    ZTCONSTSTDNIGHT4(1:NCL)=(/ 8.15,9.059,9.207,6.282,1.387,6.467,13.702,16.082,19.316,6.155,6.752,5.297,5.796,14.689,14.794,8.702,6.302,6.47,6.995,15.459,13.329,8.763,9.208,6.262,6.341,7.492,20.436  /)
    ZTCONSTAVGNIGHT5(1:NCL)=(/ -1.624,6.448,6.984,7.367,5.707,6.439,6.825,6.998,2.328,7.179,0.481,7.438,6.608,7.266,7.107,0.511,3.746,6.289,-0.175,6.379,7.448,2.427,5.872,6.396,7.122,3.064,0.139  /)
    ZTCONSTSTDNIGHT5(1:NCL)=(/ 12.051,5.187,11.018,4.543,16.577,7.822,4.511,9.389,8.91,23.547,5.818,3.847,4.502,13.879,4.403,7.281,5.518,12.255,1.407,5.101,22.598,4.451,14.495,7.386,4.51,8.287,7.079  /)
    ZTCONSTAVGNIGHT5M(1:NCL)=(/ 4.476,2.011,2.36,5.458,0.117,0.573,2.706,4.043,2.615,3.151,3.705,4.502,2.074,1.552,0.976,2.211,2.299,4.456,2.238,2.568,3.778,2.773,2.571,3.963,0.729,3.184,5.259  /)
    ZTCONSTSTDNIGHT5M(1:NCL)=(/ 12.595,1.254,4.293,2.089,7.061,4.327,6.841,2.424,5.188,6.599,4.695,11.545,1.479,5.768,4.87,3.704,2.597,2.399,1.547,9.952,8.961,5.176,1.773,15.109,6.32,10.725,1.613  /)

    ZTCONSTAVGNIGHT5G(1:NCL)=(/ 2.168,6.615,1.707,-3.289,2.13,0.324,8.895,0.584,11.117,1.814,1.293,1.714,1.647,12.19,7.416,-3.559,1.49,-4.935,1.55,5.262,1.418,4.11,1.057,2.915,1.824,1.619,0.558  /)
    ZTCONSTSTDNIGHT5G(1:NCL)=(/ 4.11,7.451,0.308,8.35,2.726,2.619,7.208,1.088,6.259,0.793,0.681,2.033,3.94,6.207,7.254,6.547,0.488,9.593,0.52,6.879,0.56,6.116,0.935,5.752,0.253,0.357,2.442  /)



! 0.3 Define constants for multi-linear regression of temp 1st layer (only ! dsn<0.20):
ZTMLRIDAY2(1:NCL)=(/2273.776,231.424,-1140.277,1077.753,-2605.152,4792.979,3877.053,393.788,&
                  &-9534.744,-11849.607,748.751,213.009,621.934,-296.89,2569.689,-870.507,&
                  &642.106,775.133,18059.695,218.558,-31549.314,-2667.837,1528.997,&
                  &-1178.686,176.003,2903.15,-2404.29 /)
ZTMLRADAY2(1:NCL)=(/-21.912,-22.981,-18.611,-3.775,19.878,-23.781,-28.105,-6.132,-39.31,&
                  &9.863,-42.609,-27.523,-7.325,-4.131,-21.878,-15.574,7.981,-6.876,5.878,&
                  &-0.976,7.58,-53.716,-14.317,-3.881,0.0,-19.631,-22.49 /)
ZTMLRBDAY2(1:NCL)=(/0.003,-0.013,0.002,0.003,0.001,0.0,-0.008,-0.001,0.003,0.001,-0.003,&
                  &0.003,0.011,-0.002,0.002,0.005,0.008,0.0,0.001,-0.0,0.001,-0.004,0.001,&
                  &0.002,-0.0,-0.004,0.001 /)
ZTMLRCDAY2(1:NCL)=(/2.63,21.656,23.678,-3.807,-0.0,-15.106,0.118,3.81,105.791,79.674,37.359,&
                  &21.693,1.89,7.147,0.852,18.914,-10.844,1.569,-136.832,0.277,225.808,73.195,&
                  &0.27,12.696,0.0,-0.986,41.659 /)
ZTMLRDDAY2(1:NCL)=(/0.046,0.043,0.043,0.007,-0.035,0.052,0.052,0.011,0.081,-0.018,0.079,0.061,&
                  &0.016,0.007,0.045,0.038,-0.019,0.013,-0.014,0.003,-0.015,0.101,0.034,&
                  &0.008,-0.0,0.036,0.041 /)
ZTMLREDAY2(1:NCL)=(/-0.002,-0.037,-0.041,0.01,0.001,0.031,0.003,-0.004,-0.192,-0.148,-0.066,&
                  &-0.037,-0.0,-0.01,0.002,-0.033,0.023,-0.0,0.254,0.0,-0.413,-0.133,0.002,&
                  &-0.021,0.001,0.004,-0.076/)

ZTMLRINIGHT2(1:NCL)=(/25.363,-4915.228,-272.045,-10385.244,-749.825,934.962,-1360.875,246.995,&
                    &-1131.312,1656.874,1945.898,-4666.717,3152.294,1125.078,15549.324,&
                    &-8619.04,-27203.797,769.829,-540.542,521.518,6036.785,127.023,&
                    &6869.104,-1328.559,885.202,-368.013,7971.363 /)
ZTMLRANIGHT2(1:NCL)=(/-21.853,-3.088,-43.653,-29.03,-16.206,1.009,4.797,6.269,-13.664,&
                    &-2.686,24.144,-11.438,-3.761,9.435,-10.072,5.866,-22.161,1.538,&
                    &-27.429,5.321,-30.082,-8.239,-7.267,-32.86,2.993,-41.127,-7.185 /)
ZTMLRBNIGHT2(1:NCL)=(/0.0,0.008,-0.007,0.002,-0.0,0.003,0.009,0.003,0.001,0.004,-0.019,&
                    &0.012,-0.001,0.002,-0.001,0.001,0.0,0.002,-0.002,-0.0,0.006,0.002,&
                    &-0.0,-0.008,-0.001,-0.004,0.0 /)
ZTMLRCNIGHT2(1:NCL)=(/21.915,40.216,45.865,103.003,22.934,-7.261,7.449,-6.425,20.528,&
                    &-9.486,-36.946,44.819,-19.352,-15.056,-103.585,59.556,220.627,-6.459,&
                    &31.715,-7.285,-18.251,7.928,-44.393,43.387,-8.759,43.923,-51.684 /)
ZTMLRDNIGHT2(1:NCL)=(/0.042,0.006,0.082,0.06,0.03,-0.003,-0.013,-0.014,0.029,0.007,-0.05,&
                    &0.026,0.007,-0.023,0.019,-0.01,0.045,-0.004,0.052,-0.013,0.064,0.015,&
                    &0.017,0.061,-0.007,0.078,0.014 /)
ZTMLRENIGHT2(1:NCL)=(/-0.039,-0.073,-0.083,-0.187,-0.041,0.017,-0.012,0.014,-0.035,0.02,0.074,&
                    &-0.081,0.039,0.03,0.192,-0.11,-0.403,0.016,-0.057,0.016,0.037,-0.012,&
                    &0.085,-0.078,0.02,-0.079,0.098/)

ZTMLRIDAY3(1:NCL)=(/-873.735,16851.732,1769.369,-1208.98,595.977,892.861,9584.742,&
                  &515.852,1914.724,-6049.516,-1895.017,-389.405,1705.265,1272.757,&
                  &29055.08,947.23,1665.529,5029.054,-333.738,24368.355,-1959.943,&
                  &1115.958,55.659,-13399.817,857.495,-1397.23,5746.104 /)
ZTMLRADAY3(1:NCL)=(/10.619,-15.266,-10.879,32.929,-25.726,14.589,-17.776,13.11,-27.152,&
                  &-42.891,14.106,-13.758,33.093,-5.929,19.702,40.24,-0.844,-20.965,&
                  &14.592,18.751,-18.933,-42.283,-39.847,7.549,21.17,-4.047,-9.189 /)
ZTMLRBDAY3(1:NCL)=(/0.0,0.0,0.0,0.001,-0.005,-0.008,-0.006,0.002,-0.004,-0.004,0.0,-0.001,&
                  &-0.003,0.004,0.002,-0.008,0.001,0.003,-0.01,-0.001,0.001,-0.008,&
                  &-0.003,0.001,-0.013,-0.001,0.0 /)
ZTMLRCDAY3(1:NCL)=(/-3.02,-109.68,-2.124,-22.655,21.564,-19.85,-53.218,-15.717,12.086,&
                  &84.146,-0.053,17.464,-44.619,-3.168,-232.189,-45.785,-11.335,-17.445,&
                  &-10.411,-196.785,32.479,33.39,39.859,90.671,-25.283,16.025,-33.291 /)
ZTMLRDDAY3(1:NCL)=(/-0.019,0.029,0.021,-0.059,0.049,-0.03,0.034,-0.027,0.052,0.087,-0.025,&
                  &0.026,-0.06,0.012,-0.039,-0.077,0.002,0.041,-0.033,-0.036,0.038,0.081,&
                  &0.075,-0.014,-0.046,0.008,0.018 /)
ZTMLREDAY3(1:NCL)=(/0.007,0.206,0.006,0.041,-0.038,0.041,0.101,0.033,-0.019,-0.153,0.002,&
                  &-0.031,0.083,0.008,0.431,0.088,0.024,0.036,0.025,0.365,-0.058,-0.06,&
                  &-0.072,-0.162,0.052,-0.029,0.064 /)

ZTMLRINIGHT3(1:27)=(/-233.194,1311.565,554.687,-170.099,685.766,5874.629,7978.278,&
                    &-8229.753,626.97,-163.629,1671.225,1383.694,22954.838,-6401.198,&
                    &1016.932,15437.012,465.208,5006.659,231532.438,-4.443,1691.159,768.856,&
                    &-1382.435,-7084.328,7818.812,1345.252,-3449.839 /)
ZTMLRANIGHT3(1:NCL)=(/-22.918,17.157,-22.351,-3.505,15.358,-58.272,-53.205,-3.579,31.338,&
                     &-5.707,15.969,23.453,28.833,-54.761,23.377,-32.093,-3.417,17.97,&
                     &11.031,-13.668,10.341,7.027,11.562,4.886,-1.373,19.383,32.732 /)
ZTMLRBNIGHT3(1:NCL)=(/-0.008,-0.002,-0.006,-0.002,-0.007,0.003,-0.008,-0.003,-0.002,&
                     &-0.007,0.009,-0.004,-0.004,-0.003,-0.006,-0.003,-0.007,-0.001,&
                     &0.017,-0.006,0.022,0.007,0.0,-0.006,-0.009,0.0,0.007 /)
ZTMLRCNIGHT3(1:NCL)=(/25.269,-25.646,16.266,5.886,-19.287,11.638,-8.116,66.062,-32.666,&
                    &7.912,-27.288,-32.353,-194.61,100.024,-29.553,-83.072,0.738,-53.752,&
                    &-1712.245,14.486,-22.205,-11.763,0.009,49.625,-55.516,-28.058,-2.126 /)
ZTMLRDNIGHT3(1:NCL)=(/0.043,-0.034,0.047,0.007,-0.029,0.115,0.104,0.006,-0.064,0.011,-0.032,&
                    &-0.045,-0.059,0.106,-0.045,0.061,0.007,-0.032,-0.021,0.026,-0.021,&
                    &-0.014,-0.021,-0.011,0.002,-0.038,-0.072 /)
ZTMLRENIGHT3(1:NCL)=(/-0.045,0.051,-0.028,-0.01,0.038,-0.019,0.018,-0.121,0.063,-0.013,&
                    &0.055,0.063,0.361,-0.182,0.058,0.158,0.0,0.099,3.15,-0.025,0.045,0.024,&
                    &0.001,-0.09,0.105,0.055,0.008 /)

! 0.4 Define centroids for Density: no distinction day and night...
! A: soT-skt B: soT-snT ; C: rsn/rsnmin
    ZRCENTRA2(1:NCL)=(/ -17.1327,-16.6176,-1406.3268,-111.5063,8582.4209,-104.4845,5945.4731,-256.0787,-123.1832,4125.603,-63.0502,116.234,93.27,&
            & -681.4148,406.5026,36.4791,272.0764,-181.7926,-40.182,8411.1836,-749.7538,-173.7347,-234.6503,150.0952,287.2344,-3390.2832,-169.1657 /)
    ZRCENTRB2(1:NCL)=(/ -227.7941,-117.5477,4.8476,-51.5537,-438.9697,-136.2922,-466.8777,159.0554,-90.4437,-726.6999,113.6042,-1473.132,-252.7198,&
            & 1629.5707,-1042.4757,54.4893,-134.9986,-19.7577,-209.336,-473.7801,-598.2844,3.2491,229.9996,-204.6059,-518.5724,-25.2679,-60.3518 /)
    ZRCENTRC2(1:NCL)=(/ 180.028,100.649,5.7877,37.0061,-508.0096,79.5346,-529.2839,-107.0202,61.9756,-825.8472,-58.4228,1093.8181,200.1835,&
            & -1316.3629,787.9306,-31.6962,-159.8319,4.9966,165.9213,-560.7889,441.1877,-6.0733,-148.6947,182.4818,406.0948,-29.5176,30.3603 /)
    ZRCENTRA3(1:NCL)=(/ 1935.7107,-339.9404,-79.7161,2301.2876,-9.5744,-642.6409,124.1115,-17680.498,-163.5648,4691.1562,-54.4623,-52.6673,-207.3275,&
            & 71.6624,-365.6533,-1692.0656,341.8825,4619.2983,-424.0497,-167.828,-26.7865,-783.2248,-69.485,-821.6069,-57.9377,308.4545,3322.5503 /)
    ZRCENTRB3(1:NCL)=(/ -190.5804,-45.5307,-250.8245,-1491.2476,-418.8027,-86.646,-434.7398,-71.8145,-286.6635,-518.4408,-119.8262,-246.0569,-172.1009,&
            & -1069.7372,-167.77,-62.9047,-856.4323,-307.5851,-630.9907,-202.6235,-248.469,-240.3025,-256.1061,-187.6713,-454.3189,-44.2783,-346.6191 /)
    ZRCENTRC3(1:NCL)=(/ -221.9388,76.8835,165.6984,1063.1085,318.6078,47.7309,259.9161,-85.7286,160.905,-603.9635,82.7302,186.4656,177.7072,&
            & 691.2839,142.9473,-76.619,556.711,-353.5813,334.5179,185.4019,119.482,74.7318,189.8985,86.954,304.084,-55.8954,-401.8718 /)
    ZRCENTRA4(1:NCL)=(/ -1146.8154,-207.3002,-180.7516,102.4712,-94.948,-711.1287,-76.8456,-281.7391,-384.2217,-300.7513,-275.5675,40.1732,-1309.3269,&
            & 1627.3419,-79.1728,-1212.6646,-309.383,-250.6138,441.7737,-308.908,-212.8091,178.9269,2907.3613,-1890.6069,1930.1302,-274.7402,-233.0804 /)
    ZRCENTRB4(1:NCL)=(/ 1220.7168,-889.175,-307.9146,-917.9799,-267.1422,-33.9146,182.684,-50.3052,-334.8376,-6.4628,-84.2195,-454.7664,890.9859,&
            & -419.7319,-922.6948,-92.3443,14.3526,139.7989,-101.8521,-112.5544,-201.776,-231.8251,-2683.9277,-46.7029,-219.9229,-125.4845,320.2653 /)
    ZRCENTRC4(1:NCL)=(/ -1041.6431,634.4794,201.7572,628.5955,201.5588,-42.6527,-91.4761,-24.9997,174.5856,-17.395,10.2183,303.4156,-711.0961,&
            & -499.0135,591.9556,-111.9868,-75.3863,-164.6775,-121.288,25.5257,108.3931,-273.3674,1781.8354,-57.3696,-271.7921,66.6956,-225.3027 /)
    ZRCENTRA5(1:NCL)=(/ -492.5684,41.3289,-518.2393,-899.2932,-1180.1403,-522.8245,-208.2451,4194.1865,-856.5204,815.509,-667.0146,-2476.4541,-544.7882,&
            & 217.9542,-292.6025,-1076.0992,174.8349,55.1778,179.227,-70.6142,-364.9809,-522.8517,607.9265,-598.3039,-1191.1541,-55.8546,-134.1052 /)
    ZRCENTRB5(1:NCL)=(/ -595.0765,-1142.4733,-501.6066,782.455,-99.0651,-630.1974,-810.9107,-329.8412,168.424,-236.2497,-77.8692,3972.958,-492.3953,&
            & -1173.1562,-687.3192,665.5585,-1064.9552,-797.7138,-266.6274,-762.4219,-796.9493,-656.3397,-206.0145,-104.9903,2015.7281,-964.4727,-867.99 /)
    ZRCENTRC5(1:NCL)=(/ 346.9941,713.2487,291.6317,-660.608,-125.1174,388.8524,493.1628,-400.1663,-214.8095,-292.1005,-101.1447,-2919.8049,284.5103,&
            & 708.0685,416.8815,-578.3391,605.3198,472.851,-330.7931,508.9077,506.0741,380.0892,-256.0961,-136.1035,-1527.0175,582.7064,523.0701 /)
    ZRCENTRA5M(1:NCL)=(/ -113.1807,-740.4465,-369.4778,-1376.5988,-183.9146,1451.9652,196.8786,900.8792,-125.511,484.7032,-247.298,750.1778,-305.5708,&
            & -1044.6616,-948.8364,-546.3726,-679.0728,110.3092,-957.2358,-1172.4614,-122.23,2076.3911,207.3752,-315.5079,33.6369,-780.0466,487.0672 /)
    ZRCENTRB5M(1:NCL)=(/ -207.8396,-44.9411,-462.5202,-101.9453,-49.9267,-1325.5876,-141.3063,-198.1673,-120.3802,-61.0082,-506.5776,-168.5689,-105.5785,&
            & -42.0027,-40.3493,-41.4848,-464.0263,-158.3677,-99.2795,-51.7693,-228.9483,-1826.0171,-202.8667,-428.9278,-1299.7432,-36.6266,-189.8887 /)
    ZRCENTRC5M(1:NCL)=(/ -5.4097,-67.1295,182.1598,-155.7934,-73.3483,705.9614,-67.0709,-245.2093,-163.3123,-88.9623,222.9218,-44.3635,-163.8482,&
            & -64.0188,-61.1647,-61.6335,208.0272,-205.9761,-150.2076,-78.9548,16.4851,1095.6368,24.8687,175.7828,659.8727,-55.3261,-238.2346 /)


    ZRCONSTAVG2(1:NCL)=(/ 0.118,0.127,0.1,0.123,0.121,0.114,0.105,0.126,0.117,0.112,0.123,0.108,0.123,&
            & 0.112,0.11,0.113,0.101,0.123,0.119,0.105,0.101,0.127,0.126,0.124,0.118,0.1,0.12 /)
    ZRCONSTSTD2(1:NCL)=(/ 0.108,0.13,0.011,0.12,0.106,0.114,0.049,0.131,0.106,0.093,0.116,0.071,0.12,&
            & 0.116,0.08,0.114,0.01,0.112,0.101,0.073,0.026,0.129,0.127,0.126,0.107,0.0,0.116 /)
    ZRCONSTAVG3(1:NCL)=(/ -2.094,-1.63,-3.117,3.317,-4.199,-1.68,-5.589,-1.587,-1.802,-1.112,-3.65,-2.939,-0.521,&
            & 1.246,-2.531,-4.27,-6.748,-2.692,-6.306,-1.527,-3.041,-9.099,-3.505,-0.899,-1.375,2.483,-1.356 /)
    ZRCONSTSTD3(1:NCL)=(/ 12.174,18.753,18.832,5.929,21.656,9.328,21.768,9.96,27.587,7.328,18.945,25.771,23.344,&
            & 13.859,14.176,13.799,22.054,14.867,17.046,25.78,22.626,17.85,21.438,10.534,24.636,21.435,8.997 /)
    ZRCONSTAVG4(1:NCL)=(/ 7.66,9.851,-3.591,-9.333,-5.518,-11.045,-6.219,-3.262,-8.867,-5.206,-4.185,-3.363,-1.635,&
            & 5.474,-9.908,-4.75,-1.167,1.248,-3.621,-1.596,-4.344,-2.691,7.249,-11.439,-2.42,-3.206,-4.395 /)
    ZRCONSTSTD4(1:NCL)=(/ 54.057,53.274,32.399,24.821,27.9,35.88,30.456,40.426,27.195,36.952,35.484,31.077,35.529,&
            & 38.284,28.679,16.02,41.395,44.715,32.746,36.108,34.602,30.86,43.14,29.05,37.5,33.95,33.144 /)
    ZRCONSTAVG5(1:NCL)=(/ 7.355,4.431,7.464,1.963,-12.431,10.64,10.003,2.731,2.882,5.387,-7.979,2.501,8.535,&
            & 0.361,7.241,3.097,20.792,6.269,0.698,10.866,13.567,6.16,0.221,-8.132,1.602,-1.968,9.121 /)
    ZRCONSTSTD5(1:NCL)=(/ 31.357,38.082,34.021,28.494,17.254,34.925,36.704,28.34,25.373,38.297,18.477,30.197,33.887,&
            & 35.295,33.583,27.921,46.365,36.508,34.124,35.693,39.794,29.518,26.088,23.593,27.908,29.329,37.779 /)
    ZRCONSTAVG5M(1:NCL)=(/ 1.831,2.629,1.608,2.097,1.145,5.817,1.771,6.172,4.617,2.478,1.753,2.836,3.415,&
            & 2.396,4.043,1.64,2.303,6.267,0.652,1.863,1.305,4.164,2.017,2.307,7.068,2.227,5.421 /)
    ZRCONSTSTD5M(1:NCL)=(/ 16.309,15.137,14.247,20.153,14.305,22.913,12.115,27.239,21.127,15.187,15.076,14.407,28.652,&
            & 17.864,19.886,13.371,15.589,25.989,19.55,18.851,12.228,14.894,14.335,12.719,22.919,15.947,27.081 /)


    ZRMLRI2(1:NCL)=(/ -17.1327,-16.6176,-1406.3268,-111.5063,8582.4209,-104.4845,5945.4731,-256.0787,-123.1832,4125.603,-63.0502,116.234,93.27,&
            & -681.4148,406.5026,36.4791,272.0764,-181.7926,-40.182,8411.1836,-749.7538,-173.7347,-234.6503,150.0952,287.2344,-3390.2832,-169.1657 /)
    ZRMLRA2(1:NCL)=(/ -227.7941,-117.5477,4.8476,-51.5537,-438.9697,-136.2922,-466.8777,159.0554,-90.4437,-726.6999,113.6042,-1473.132,-252.7198,&
            & 1629.5707,-1042.4757,54.4893,-134.9986,-19.7577,-209.336,-473.7801,-598.2844,3.2491,229.9996,-204.6059,-518.5724,-25.2679,-60.3518 /)
    ZRMLRB2(1:NCL)=(/ 0.9961,2.4539,0.6722,1.1979,0.365,2.1321,3.4408,1.6125,1.1548,2.276,2.0354,1.0142,0.7555,&
            & 2.0656,2.3182,2.6313,0.1691,1.1362,1.1895,-2.2311,0.8056,1.1813,2.274,2.5718,1.3692,0.8091,0.7002 /)
    ZRMLRC2(1:NCL)=(/ 0.2917,-0.3529,5.2271,0.3738,-29.5811,0.1503,-21.8666,0.4857,0.4921,-13.2726,-0.3882,1.3774,-0.0067,&
            & 0.2359,-0.7273,-0.8857,-0.0009,0.6388,0.2861,-27.2512,3.5548,0.5669,0.0649,-0.9862,-0.6329,12.5992,0.8081 /)
    ZRMLRD2(1:NCL)=(/ 0.0003,-0.0038,0.0012,-0.0001,0.0004,-0.0025,-0.0034,-0.0011,-0.0001,-0.0025,-0.0022,0.0,0.0009,&
            & -0.0019,-0.0031,-0.004,0.0018,-0.0,-0.0001,0.004,0.0007,-0.0002,-0.0029,-0.004,-0.0006,0.0008,0.001 /)
    ZRMLRE2(1:NCL)=(/ 180.028,100.649,5.7877,37.0061,-508.0096,79.5346,-529.2839,-107.0202,61.9756,-825.8472,-58.4228,1093.8181,200.1835,&
            & -1316.3629,787.9306,-31.6962,-159.8319,4.9966,165.9213,-560.7889,441.1877,-6.0733,-148.6947,182.4818,406.0948,-29.5176,30.3603 /)
    ZRMLRI3(1:NCL)=(/ 1935.7107,-339.9404,-79.7161,2301.2876,-9.5744,-642.6409,124.1115,-17680.498,-163.5648,4691.1562,-54.4623,-52.6673,-207.3275,&
            & 71.6624,-365.6533,-1692.0656,341.8825,4619.2983,-424.0497,-167.828,-26.7865,-783.2248,-69.485,-821.6069,-57.9377,308.4545,3322.5503 /)
    ZRMLRA3(1:NCL)=(/ -190.5804,-45.5307,-250.8245,-1491.2476,-418.8027,-86.646,-434.7398,-71.8145,-286.6635,-518.4408,-119.8262,-246.0569,-172.1009,&
            & -1069.7372,-167.77,-62.9047,-856.4323,-307.5851,-630.9907,-202.6235,-248.469,-240.3025,-256.1061,-187.6713,-454.3189,-44.2783,-346.6191 /)
    ZRMLRB3(1:NCL)=(/ 2.5887,3.4935,-0.2986,-26.9898,1.0548,2.5842,-1.1878,-0.2668,-0.3841,6.0931,0.3392,1.0216,2.5854,&
            & -2.4866,3.2367,0.7124,-1.0935,-0.0613,-1.8165,3.2643,-0.7194,-1.4905,0.4214,2.2662,-0.412,0.5233,1.5688 /)
    ZRMLRC3(1:NCL)=(/ -7.4723,0.3396,1.12,0.2602,0.4694,2.0554,0.9941,65.5369,1.5684,-18.8703,0.5914,0.4567,0.2375,&
            & 2.1997,0.7521,6.5562,0.7408,-15.3414,3.6793,0.0003,1.1461,4.3173,0.762,3.1453,1.3497,-0.7248,-11.3266 /)
    ZRMLRD3(1:NCL)=(/ -0.0022,-0.0063,0.0029,0.1035,0.0002,-0.004,0.0049,0.0029,0.003,-0.007,0.0015,-0.0,-0.0033,&
            & 0.0133,-0.0057,0.0007,0.0044,0.0015,0.0059,-0.006,0.0038,0.0053,0.0015,-0.0033,0.0032,0.0009,-0.001 /)
    ZRMLRE3(1:NCL)=(/ -221.9388,76.8835,165.6984,1063.1085,318.6078,47.7309,259.9161,-85.7286,160.905,-603.9635,82.7302,186.4656,177.7072,&
            & 691.2839,142.9473,-76.619,556.711,-353.5813,334.5179,185.4019,119.482,74.7318,189.8985,86.954,304.084,-55.8954,-401.8718 /)
    ZRMLRI4(1:NCL)=(/ -1146.8154,-207.3002,-180.7516,102.4712,-94.948,-711.1287,-76.8456,-281.7391,-384.2217,-300.7513,-275.5675,40.1732,-1309.3269,&
            & 1627.3419,-79.1728,-1212.6646,-309.383,-250.6138,441.7737,-308.908,-212.8091,178.9269,2907.3613,-1890.6069,1930.1302,-274.7402,-233.0804 /)
    ZRMLRA4(1:NCL)=(/ 1220.7168,-889.175,-307.9146,-917.9799,-267.1422,-33.9146,182.684,-50.3052,-334.8376,-6.4628,-84.2195,-454.7664,890.9859,&
            & -419.7319,-922.6948,-92.3443,14.3526,139.7989,-101.8521,-112.5544,-201.776,-231.8251,-2683.9277,-46.7029,-219.9229,-125.4845,320.2653 /)
    ZRMLRB4(1:NCL)=(/ 3.2574,0.6466,1.8436,-0.7315,3.474,-1.7113,3.9412,-0.5139,-0.8465,-0.1288,-0.1337,2.3852,4.7711,&
            & 0.7581,-2.0799,-0.4427,-0.5633,-1.3495,0.6653,2.3494,0.4349,0.9131,-28.8256,0.4061,-0.3226,0.7883,3.9233 /)
    ZRMLRC4(1:NCL)=(/ 2.3842,2.0207,0.79,1.536,-0.2226,3.8077,-1.1879,1.8054,2.6854,1.5804,1.6633,-0.0656,2.8824,&
            & -4.2799,2.7807,5.4022,1.8703,1.774,-1.0592,0.9242,1.3474,0.1538,0.2848,7.3374,-5.581,1.3074,-0.6656 /)
    ZRMLRD4(1:NCL)=(/ -0.006,0.0007,-0.0024,0.0032,-0.0068,0.0059,-0.0076,0.003,0.0037,0.0026,0.0023,-0.0038,-0.0113,&
            & 0.0,0.0065,0.0031,0.0031,0.005,0.0004,-0.0039,0.0008,0.0001,0.1099,0.0013,0.0019,0.0002,-0.0078 /)
    ZRMLRE4(1:NCL)=(/ -1041.6431,634.4794,201.7572,628.5955,201.5588,-42.6527,-91.4761,-24.9997,174.5856,-17.395,10.2183,303.4156,-711.0961,&
            & -499.0135,591.9556,-111.9868,-75.3863,-164.6775,-121.288,25.5257,108.3931,-273.3674,1781.8354,-57.3696,-271.7921,66.6956,-225.3027 /)
    ZRMLRI5(1:NCL)=(/ -492.5684,41.3289,-518.2393,-899.2932,-1180.1403,-522.8245,-208.2451,4194.1865,-856.5204,815.509,-667.0146,-2476.4541,-544.7882,&
            & 217.9542,-292.6025,-1076.0992,174.8349,55.1778,179.227,-70.6142,-364.9809,-522.8517,607.9265,-598.3039,-1191.1541,-55.8546,-134.1052 /)
    ZRMLRA5(1:NCL)=(/ -595.0765,-1142.4733,-501.6066,782.455,-99.0651,-630.1974,-810.9107,-329.8412,168.424,-236.2497,-77.8692,3972.958,-492.3953,&
            & -1173.1562,-687.3192,665.5585,-1064.9552,-797.7138,-266.6274,-762.4219,-796.9493,-656.3397,-206.0145,-104.9903,2015.7281,-964.4727,-867.99 /)
    ZRMLRB5(1:NCL)=(/ 3.4988,-0.6269,2.205,4.0731,0.3353,2.2156,-0.6265,0.9001,4.4237,-1.7768,-0.1871,8.1571,3.3276,&
            & -1.8436,1.0559,4.5966,-0.7507,-1.2074,-1.0449,2.186,1.1034,3.7849,-0.9617,-1.5213,3.9069,-0.9368,-0.194 /)
    ZRMLRC5(1:NCL)=(/ 2.0361,2.2459,2.3957,1.7521,4.964,2.5651,2.7524,-14.4073,2.2163,-0.5034,3.2953,2.5924,2.1102,&
            & 2.1178,2.2227,2.4199,1.8494,1.9398,1.5512,0.9832,2.5829,2.1705,-0.3938,3.7468,1.2575,2.5191,2.3563 /)
    ZRMLRD5(1:NCL)=(/ -0.0079,0.0023,-0.0045,-0.0105,0.0013,-0.0043,0.0023,0.0004,-0.011,0.0039,0.002,-0.0254,-0.0074,&
            & 0.0052,-0.0016,-0.0118,0.0024,0.0038,0.0028,-0.0044,-0.0014,-0.0087,0.0029,0.0047,-0.0104,0.003,0.0015 /)
    ZRMLRE5(1:NCL)=(/ 346.9941,713.2487,291.6317,-660.608,-125.1174,388.8524,493.1628,-400.1663,-214.8095,-292.1005,-101.1447,-2919.8049,284.5103,&
            & 708.0685,416.8815,-578.3391,605.3198,472.851,-330.7931,508.9077,506.0741,380.0892,-256.0961,-136.1035,-1527.0175,582.7064,523.0701 /)
    ZRMLRI5M(1:NCL)=(/ -113.1807,-740.4465,-369.4778,-1376.5988,-183.9146,1451.9652,196.8786,900.8792,-125.511,484.7032,-247.298,750.1778,-305.5708,&
            & -1044.6616,-948.8364,-546.3726,-679.0728,110.3092,-957.2358,-1172.4614,-122.23,2076.3911,207.3752,-315.5079,33.6369,-780.0466,487.0672 /)
    ZRMLRA5M(1:NCL)=(/ -207.8396,-44.9411,-462.5202,-101.9453,-49.9267,-1325.5876,-141.3063,-198.1673,-120.3802,-61.0082,-506.5776,-168.5689,-105.5785,&
            & -42.0027,-40.3493,-41.4848,-464.0263,-158.3677,-99.2795,-51.7693,-228.9483,-1826.0171,-202.8667,-428.9278,-1299.7432,-36.6266,-189.8887 /)
    ZRMLRB5M(1:NCL)=(/ 0.0266,0.7795,0.4634,1.4432,-0.2284,-2.7732,-0.4469,-2.1434,-2.163,-0.4412,0.1634,-0.6046,0.7437,&
            & 1.7888,1.3795,0.5039,0.6668,-2.3508,1.839,2.0634,-0.2569,-2.6629,-0.6063,-0.0142,-2.5617,1.3443,-2.2609 /)
    ZRMLRC5M(1:NCL)=(/ 1.5099,3.3012,2.6351,5.9815,1.5899,-1.1878,0.5619,-0.8123,2.6659,-0.7224,2.3194,-1.4155,2.0163,&
            & 4.0086,3.8651,2.6144,3.6931,2.091,4.1759,4.5137,1.728,-2.9971,0.5366,2.5993,3.9186,3.1262,0.8369 /)
    ZRMLRD5M(1:NCL)=(/ 0.0005,-0.0012,-0.0005,-0.0027,0.0009,0.0053,0.0012,0.0046,0.0045,0.0014,0.0001,0.0016,-0.0005,&
            & -0.0033,-0.0026,-0.0006,-0.001,0.0048,-0.0036,-0.004,0.0009,0.0051,0.0016,0.0004,0.005,-0.0024,0.0046 /)
    ZRMLRE5M(1:NCL)=(/ -5.4097,-67.1295,182.1598,-155.7934,-73.3483,705.9614,-67.0709,-245.2093,-163.3123,-88.9623,222.9218,-44.3635,-163.8482,&
            & -64.0188,-61.1647,-61.6335,208.0272,-205.9761,-150.2076,-78.9548,16.4851,1095.6368,24.8687,175.7828,659.8727,-55.3261,-238.2346 /)


!*******************************************************************************
! 2. Start snow parametrizations
!    Here we start the snow temperature and density parametrisation:
!    both are simple exponential function. Temperature relaxes to first soil layer
!    at the bottom, while the density relaxes to mean density. 
!    The missing snow mass it is added to the bottom layer simply increasing the
!    density of this layer keeping the depth fixed.
!*******************************************************************************
  DO JL=KIDIA, KFDIA
    IF (PMU0(JL) > ZEPSILON ) THEN ! Daytime
        ! Initialise values to avoid floating point errors
        ZTCENTRA=ZTCENTRADAY2
        ZTCENTRB=ZTCENTRBDAY2
        ZTCENTRC=ZTCENTRCDAY2
        ZTCONSTAVG=ZTCONSTAVGDAY2
        ZTCONSTSTD=ZTCONSTSTDDAY2
        KLEVMID(JL)=MAX(KLEVSNA(JL)-1,1)

    IF ( PSSN(JL) < ZSNPERT .AND. LDLAND(JL) ) THEN ! seasonal snow
      IF ( ZDSNTOT(JL) < 0.15_JPRB ) THEN 
        ZTCENTRA=ZTCENTRADAY2
        ZTCENTRB=ZTCENTRBDAY2
        ZTCENTRC=ZTCENTRCDAY2

        ZTCONSTAVG=ZTCONSTAVGDAY2
        ZTCONSTSTD=ZTCONSTSTDDAY2

        ZTMLRI=ZTMLRIDAY2
        ZTMLRA=ZTMLRADAY2
        ZTMLRB=ZTMLRBDAY2
        ZTMLRC=ZTMLRCDAY2
        ZTMLRD=ZTMLRDDAY2
        ZTMLRE=ZTMLREDAY2

      ELSEIF ( ZDSNTOT(JL) < 0.20_JPRB ) THEN 
        ZTCENTRA=ZTCENTRADAY3
        ZTCENTRB=ZTCENTRBDAY3
        ZTCENTRC=ZTCENTRCDAY3

        ZTCONSTAVG=ZTCONSTAVGDAY3
        ZTCONSTSTD=ZTCONSTSTDDAY3
        
        ZTMLRI=ZTMLRIDAY3
        ZTMLRA=ZTMLRADAY3
        ZTMLRB=ZTMLRBDAY3
        ZTMLRC=ZTMLRCDAY3
        ZTMLRD=ZTMLRDDAY3
        ZTMLRE=ZTMLREDAY3

      ELSEIF ( ZDSNTOT(JL) < 0.25_JPRB ) THEN 
        ZTCENTRA=ZTCENTRADAY4
        ZTCENTRB=ZTCENTRBDAY4
        ZTCENTRC=ZTCENTRCDAY4
        
        ZTCONSTAVG=ZTCONSTAVGDAY4
        ZTCONSTSTD=ZTCONSTSTDDAY4

      ELSEIF ( ZDSNTOT(JL) < 0.50_JPRB ) THEN 
        IF (PSDOR(JL)<ZSDORTHR)THEN
          ZTCENTRA=ZTCENTRADAY5
          ZTCENTRB=ZTCENTRBDAY5
          ZTCENTRC=ZTCENTRCDAY5

          ZTCONSTAVG=ZTCONSTAVGDAY5
          ZTCONSTSTD=ZTCONSTSTDDAY5
        ELSEIF (KLEVSNA(JL)==2)THEN
          ZTCENTRA=ZTCENTRADAY2
          ZTCENTRB=ZTCENTRBDAY2
          ZTCENTRC=ZTCENTRCDAY2

          ZTCONSTAVG=ZTCONSTAVGDAY2
          ZTCONSTSTD=ZTCONSTSTDDAY2
        ELSEIF (KLEVSNA(JL)==3)THEN
          ZTCENTRA=ZTCENTRADAY3
          ZTCENTRB=ZTCENTRBDAY3
          ZTCENTRC=ZTCENTRCDAY3

          ZTCONSTAVG=ZTCONSTAVGDAY3
          ZTCONSTSTD=ZTCONSTSTDDAY3

        ELSEIF (KLEVSNA(JL)==4)THEN
          ZTCENTRA=ZTCENTRADAY4
          ZTCENTRB=ZTCENTRBDAY4
          ZTCENTRC=ZTCENTRCDAY4

          ZTCONSTAVG=ZTCONSTAVGDAY4
          ZTCONSTSTD=ZTCONSTSTDDAY4
        ELSEIF (KLEVSNA(JL)==5)THEN
          ZTCENTRA=ZTCENTRADAY5
          ZTCENTRB=ZTCENTRBDAY5
          ZTCENTRC=ZTCENTRCDAY5

          ZTCONSTAVG=ZTCONSTAVGDAY5
          ZTCONSTSTD=ZTCONSTSTDDAY5
        ENDIF
      ELSEIF ( ZDSNTOT(JL) >= 0.50_JPRB ) THEN 
        ZTCENTRA=ZTCENTRADAY5M
        ZTCENTRB=ZTCENTRBDAY5M
        ZTCENTRC=ZTCENTRCDAY5M

        ZTCONSTAVG=ZTCONSTAVGDAY5M
        ZTCONSTSTD=ZTCONSTSTDDAY5M
      ENDIF
    ELSE ! Glaciers, daytime
        ZTCENTRA=ZTCENTRADAY5G
        ZTCENTRB=ZTCENTRBDAY5G
        ZTCENTRC=ZTCENTRCDAY5G

        ZTCONSTAVG=ZTCONSTAVGDAY5G
        ZTCONSTSTD=ZTCONSTSTDDAY5G
    ENDIF
    ELSEIF (PMU0(JL)<=ZEPSILON) THEN !nighttime
        ! Initialise values to avoid floating point errors
        ZTCENTRA=ZTCENTRANIGHT2
        ZTCENTRB=ZTCENTRBNIGHT2
        ZTCENTRC=ZTCENTRCNIGHT2
        ZTCONSTAVG=ZTCONSTAVGNIGHT2
        ZTCONSTSTD=ZTCONSTSTDNIGHT2
        KLEVMID(JL)=MAX(KLEVSNA(JL)-1,1)

    IF ( PSSN(JL) < ZSNPERT .AND. LDLAND(JL) ) THEN ! seasonal snow
      IF ( ZDSNTOT(JL) < 0.15_JPRB ) THEN 
        ZTCENTRA=ZTCENTRANIGHT2
        ZTCENTRB=ZTCENTRBNIGHT2
        ZTCENTRC=ZTCENTRCNIGHT2

        ZTCONSTAVG=ZTCONSTAVGNIGHT2
        ZTCONSTSTD=ZTCONSTSTDNIGHT2

        ZTMLRI=ZTMLRINIGHT2
        ZTMLRA=ZTMLRANIGHT2
        ZTMLRB=ZTMLRBNIGHT2
        ZTMLRC=ZTMLRCNIGHT2
        ZTMLRD=ZTMLRDNIGHT2
        ZTMLRE=ZTMLRENIGHT2
        
      ELSEIF ( ZDSNTOT(JL) < 0.20_JPRB ) THEN 
        ZTCENTRA=ZTCENTRANIGHT3
        ZTCENTRB=ZTCENTRBNIGHT3
        ZTCENTRC=ZTCENTRCNIGHT3

        ZTCONSTAVG=ZTCONSTAVGNIGHT3
        ZTCONSTSTD=ZTCONSTSTDNIGHT3
        
        ZTMLRI=ZTMLRINIGHT3
        ZTMLRA=ZTMLRANIGHT3
        ZTMLRB=ZTMLRBNIGHT3
        ZTMLRC=ZTMLRCNIGHT3
        ZTMLRD=ZTMLRDNIGHT3
        ZTMLRE=ZTMLRENIGHT3

      ELSEIF ( ZDSNTOT(JL) < 0.25_JPRB ) THEN 
        ZTCENTRA=ZTCENTRANIGHT4
        ZTCENTRB=ZTCENTRBNIGHT4
        ZTCENTRC=ZTCENTRCNIGHT4

        ZTCONSTAVG=ZTCONSTAVGNIGHT4
        ZTCONSTSTD=ZTCONSTSTDNIGHT4
        
      ELSEIF ( ZDSNTOT(JL) < 0.50_JPRB ) THEN 
        IF (PSDOR(JL)<ZSDORTHR)THEN
          ZTCENTRA=ZTCENTRANIGHT5
          ZTCENTRB=ZTCENTRBNIGHT5
          ZTCENTRC=ZTCENTRCNIGHT5

          ZTCONSTAVG=ZTCONSTAVGNIGHT5
          ZTCONSTSTD=ZTCONSTSTDNIGHT5
        ELSEIF (KLEVSNA(JL)==2)THEN
          ZTCENTRA=ZTCENTRANIGHT2
          ZTCENTRB=ZTCENTRBNIGHT2
          ZTCENTRC=ZTCENTRCNIGHT2

          ZTCONSTAVG=ZTCONSTAVGNIGHT2
          ZTCONSTSTD=ZTCONSTSTDNIGHT2
        ELSEIF (KLEVSNA(JL)==3)THEN
          ZTCENTRA=ZTCENTRANIGHT3
          ZTCENTRB=ZTCENTRBNIGHT3
          ZTCENTRC=ZTCENTRCNIGHT3

          ZTCONSTAVG=ZTCONSTAVGNIGHT3
          ZTCONSTSTD=ZTCONSTSTDNIGHT3
        ELSEIF (KLEVSNA(JL)==4)THEN
          ZTCENTRA=ZTCENTRANIGHT4
          ZTCENTRB=ZTCENTRBNIGHT4
          ZTCENTRC=ZTCENTRCNIGHT4

          ZTCONSTAVG=ZTCONSTAVGNIGHT4
          ZTCONSTSTD=ZTCONSTSTDNIGHT4
        ELSEIF (KLEVSNA(JL)==5)THEN
          ZTCENTRA=ZTCENTRANIGHT5
          ZTCENTRB=ZTCENTRBNIGHT5
          ZTCENTRC=ZTCENTRCNIGHT5

          ZTCONSTAVG=ZTCONSTAVGNIGHT5
          ZTCONSTSTD=ZTCONSTSTDNIGHT5
        ENDIF
        
      ELSEIF ( ZDSNTOT(JL) >= 0.50_JPRB ) THEN 
        ZTCENTRA=ZTCENTRANIGHT5M
        ZTCENTRB=ZTCENTRBNIGHT5M
        ZTCENTRC=ZTCENTRCNIGHT5M

        ZTCONSTAVG=ZTCONSTAVGNIGHT5M
        ZTCONSTSTD=ZTCONSTSTDNIGHT5M
      ENDIF
    ELSE !Glaciers, nighttime
        ZTCENTRA=ZTCENTRANIGHT5G
        ZTCENTRB=ZTCENTRBNIGHT5G
        ZTCENTRC=ZTCENTRCNIGHT5G

        ZTCONSTAVG=ZTCONSTAVGNIGHT5G
        ZTCONSTSTD=ZTCONSTSTDNIGHT5G
    ENDIF
    ENDIF

    PTCONSTAVG(JL,1:NCL)=ZTCONSTAVG(1:NCL)
    PTCONSTSTD(JL,1:NCL)=ZTCONSTSTD(1:NCL)

! Assign density depending on snow depth::
    ! Initialise arrays, to avoid undesired effects.
    ZRMLRA(1:NCL)=ZRMLRA2(1:NCL)
    ZRMLRB(1:NCL)=ZRMLRB2(1:NCL)
    ZRMLRC(1:NCL)=ZRMLRC2(1:NCL)
    ZRMLRD(1:NCL)=ZRMLRD2(1:NCL)
    ZRMLRE(1:NCL)=ZRMLRE2(1:NCL)
    IF ( ZDSNTOT(JL) < 0.15_JPRB ) THEN 
      ZRCENTRA=ZRCENTRA2
      ZRCENTRB=ZRCENTRB2
      ZRCENTRC=ZRCENTRC2

      ZRCONSTAVG=ZRCONSTAVG2
      ZRCONSTSTD=ZRCONSTSTD2

      ZRMLRI=ZRMLRI2
      ZRMLRA=ZRMLRA2
      ZRMLRB=ZRMLRB2
      ZRMLRC=ZRMLRC2
      ZRMLRD=ZRMLRD2
      ZRMLRE=ZRMLRE2
      
      KLEVMID(JL)=2_JPIM
    ELSEIF ( ZDSNTOT(JL) < 0.20_JPRB ) THEN 
      ZRCENTRA=ZRCENTRA3
      ZRCENTRB=ZRCENTRB3
      ZRCENTRC=ZRCENTRC3

      ZRCONSTAVG=ZRCONSTAVG3
      ZRCONSTSTD=ZRCONSTSTD3
      
      ZRMLRI=ZRMLRI3
      ZRMLRA=ZRMLRA3
      ZRMLRB=ZRMLRB3
      ZRMLRC=ZRMLRC3
      ZRMLRD=ZRMLRD3
      ZRMLRE=ZRMLRE3

      KLEVMID(JL)=2_JPIM

    ELSEIF ( ZDSNTOT(JL) < 0.25_JPRB ) THEN 
      ZRCENTRA=ZRCENTRA4
      ZRCENTRB=ZRCENTRB4
      ZRCENTRC=ZRCENTRC4

      ZRCONSTAVG=ZRCONSTAVG4
      ZRCONSTSTD=ZRCONSTSTD4
      
      ZRMLRI=ZRMLRI4
      ZRMLRA=ZRMLRA4
      ZRMLRB=ZRMLRB4
      ZRMLRC=ZRMLRC4
      ZRMLRD=ZRMLRD4
      ZRMLRE=ZRMLRE4
 
      KLEVMID(JL)=3_JPIM

    ELSEIF ( ZDSNTOT(JL) < 0.50_JPRB ) THEN 
      IF (PSDOR(JL)<ZSDORTHR)THEN
        ZRCENTRA=ZRCENTRA5
        ZRCENTRB=ZRCENTRB5
        ZRCENTRC=ZRCENTRC5

        ZRCONSTAVG=ZRCONSTAVG5
        ZRCONSTSTD=ZRCONSTSTD5
        
        ZRMLRI=ZRMLRI5
        ZRMLRA=ZRMLRA5
        ZRMLRB=ZRMLRB5
        ZRMLRC=ZRMLRC5
        ZRMLRD=ZRMLRD5
        ZRMLRE=ZRMLRE5
      ELSEIF (KLEVSNA(JL)==2)THEN
        ZRCENTRA=ZRCENTRA2
        ZRCENTRB=ZRCENTRB2
        ZRCENTRC=ZRCENTRC2

        ZRCONSTAVG=ZRCONSTAVG2
        ZRCONSTSTD=ZRCONSTSTD2
        
        ZRMLRI=ZRMLRI2
        ZRMLRA=ZRMLRA2
        ZRMLRB=ZRMLRB2
        ZRMLRC=ZRMLRC2
        ZRMLRD=ZRMLRD2
        ZRMLRE=ZRMLRE2
      ELSEIF (KLEVSNA(JL)==3)THEN
        ZRCENTRA=ZRCENTRA3
        ZRCENTRB=ZRCENTRB3
        ZRCENTRC=ZRCENTRC3

        ZRCONSTAVG=ZRCONSTAVG3
        ZRCONSTSTD=ZRCONSTSTD3
        
        ZRMLRI=ZRMLRI3
        ZRMLRA=ZRMLRA3
        ZRMLRB=ZRMLRB3
        ZRMLRC=ZRMLRC3
        ZRMLRD=ZRMLRD3
        ZRMLRE=ZRMLRE3
      ELSEIF (KLEVSNA(JL)==4)THEN
        ZRCENTRA=ZRCENTRA4
        ZRCENTRB=ZRCENTRB4
        ZRCENTRC=ZRCENTRC4

        ZRCONSTAVG=ZRCONSTAVG4
        ZRCONSTSTD=ZRCONSTSTD4
        
        ZRMLRI=ZRMLRI4
        ZRMLRA=ZRMLRA4
        ZRMLRB=ZRMLRB4
        ZRMLRC=ZRMLRC4
        ZRMLRD=ZRMLRD4
        ZRMLRE=ZRMLRE4
      ELSEIF (KLEVSNA(JL)==5)THEN
        ZRCENTRA=ZRCENTRA5
        ZRCENTRB=ZRCENTRB5
        ZRCENTRC=ZRCENTRC5

        ZRCONSTAVG=ZRCONSTAVG5
        ZRCONSTSTD=ZRCONSTSTD5
        
        ZRMLRI=ZRMLRI5
        ZRMLRA=ZRMLRA5
        ZRMLRB=ZRMLRB5
        ZRMLRC=ZRMLRC5
        ZRMLRD=ZRMLRD5
        ZRMLRE=ZRMLRE5
      ENDIF

      IF (PSDOR(JL)<ZSDORTHR)THEN
        KLEVMID(JL)=3_JPIM
      ELSEIF (KLEVSNA(JL)<KLEVSN .AND. KLEVSNA(JL)>2)THEN
        KLEVMID(JL)=KLEVSNA(JL)-1
      ELSEIF (KLEVSNA(JL)==2)THEN
        KLEVMID(JL)=KLEVSNA(JL)
      ELSE
        KLEVMID(JL)=3_JPIM
      ENDIF

    ELSEIF ( ZDSNTOT(JL) >= 0.50_JPRB ) THEN 
      ZRCENTRA=ZRCENTRA5M
      ZRCENTRB=ZRCENTRB5M
      ZRCENTRC=ZRCENTRC5M

      ZRCONSTAVG=ZRCONSTAVG5M
      ZRCONSTSTD=ZRCONSTSTD5M
      
      ZRMLRI=ZRMLRI5M
      ZRMLRA=ZRMLRA5M
      ZRMLRB=ZRMLRB5M
      ZRMLRC=ZRMLRC5M
      ZRMLRD=ZRMLRD5M
      ZRMLRE=ZRMLRE5M

      KLEVMID(JL)=MAX(KLEVSNA(JL)-1,1)
     
    ENDIF
! Final assignment:
!  RCENTRA(JL)=ZRCENTRA
!  RCENTRB(JL)=ZRCENTRB
!  RCENTRC(JL)=ZRCENTRC

   PRCONSTAVG(JL,1:NCL)=ZRCONSTAVG(1:NCL)
   PRCONSTSTD(JL,1:NCL)=ZRCONSTSTD(1:NCL)
   
!  RMLRI(JL)=ZRMLRI
!  RMLRA(JL)=ZRMLRA
!  RMLRB(JL)=ZRMLRB
!  RMLRC(JL)=ZRMLRC
!  RMLRD(JL)=ZRMLRD
!  RMLRE(JL)=ZRMLRE

!*****************************************
! 2.1 Find closest cluster centre using euclidean metric:
!     Features are common to temp and dens
   ZFEATA=PTSOIL(JL)-PTSKIN(JL)
   ZFEATB=PTSOIL(JL)-PTSN(JL)
   ZFEATC=PRSN(JL)/RHOMINSND

! 2.1.1 Temperature
   DO KCL=1,NCL
     ZTDIST(KCL)=SQRT( (ZFEATA-ZTCENTRA(KCL))**2_JPRB + (ZFEATB-ZTCENTRB(KCL))**2_JPRB + (ZFEATC-ZTCENTRC(KCL))**2_JPRB )
   ENDDO
   PTMINCL(JL)=MINLOC(ZTDIST,DIM=1)

! 2.1.2 Density:
   DO KCL=1,NCL
     ZRDIST(KCL)=SQRT( (ZFEATA-ZRCENTRA(KCL))**2_JPRB + (ZFEATB-ZRCENTRB(KCL))**2_JPRB + (ZFEATC-ZRCENTRC(KCL))**2_JPRB )
   ENDDO
   PRMINCL(JL)=MINLOC(ZRDIST,DIM=1)

!*****************************************
! 2.1 Initialize warm start (WS) variables
    IF ( PSSN(JL) < ZSNPERT .AND. LDLAND(JL)) THEN
        PTSNWS(JL, 1)          = MIN(RTT,PTSKIN(JL))
        PTSNWS(JL, 2:KLEVSN-1) = PTSN(JL)
        PTSNWS(JL, KLEVSN)     = MIN(RTT,PTSOIL(JL))

        PRSNWS(JL, 1:KLEVSN) = PRSN(JL)
        PRSNMAX(JL)          = RHOMAXSN_NEW

        PSSNWS(JL, 1)        = PSSN(JL)
        PWSNWS(JL, 1)        = PWSN(JL)
        PSSNWS(JL, 2:KLEVSN) = 0._JPRB
        PWSNWS(JL, 2:KLEVSN) = 0._JPRB

!----- Density
        ZP1=ZRMLRI(PRMINCL(JL))+ZRMLRA(PRMINCL(JL))*PASN(JL)
        ZP2=                    ZRMLRB(PRMINCL(JL))*PRSN(JL)
        ZP3=                    ZRMLRC(PRMINCL(JL))*PTSN(JL)
        ZP4=                    ZRMLRD(PRMINCL(JL))*(PRSN(JL))**2_JPRB
        ZP5=                    ZRMLRE(PRMINCL(JL))*(PASN(JL))**2_JPRB

        PRSNTOP(JL)=ZP1+ZP2+ZP3+ZP4+ZP5

        PTSNBOTTOM(JL)=MIN(RTT,PTSOIL(JL))
        PTSNMIDDLE(JL)=MIN(RTT,0.5_JPRB*(PTSN(JL)+PTSOIL(JL)))
        IF ( ZDSNTOT(JL) < 0.20_JPRB ) THEN 
          PTSNTOP(JL) = PTSN(JL)
        ELSE
          PTSNTOP(JL) = PTSKIN(JL)
        ENDIF
        IF ( ZDSNTOT(JL) < 0.20_JPRB ) THEN
          ZSADEPTH(JL)=0._JPRB
        ELSE
          ZSADEPTH(JL)=0.5_JPRB*RDAT(1)
          IF (PTSNTOP(JL)>=PTSN(JL)) THEN
            ZACTDEPTH(JL)=ZSNDEPTH(JL, KLEVSNA(JL)-1)
          ELSE
            ZACTDEPTH(JL)=ZDSNTOT(JL)
          ENDIF
        ENDIF

    ELSE ! Glacier or sea-ice ini
        KLEVSNA(JL)=KLEVSN !KSNACC+1

        PTSNWS(JL, 1)        = MIN(RTT,PTSKIN(JL))
        PTSNWS(JL, 2:KLEVSN) = MIN(RTT,PTSN(JL))

        PRSNMAX(JL)              = 300._JPRB
        PRSNWS(JL,1:KLEVSN)      = PRSNMAX(JL)
       !PRSNWS(JL,KSNACC:KLEVSN) = PRSNMAX(JL)
        PSSNWS(JL,1:KLEVSN)      = PRSNMAX(JL)*ZDSNR(JL,1:KLEVSN)
        PWSNWS(JL,1:KLEVSN)      = 0._JPRB
!----- Active depth
      PTSNTOP(JL)   = PTSKIN(JL)
      IF (PTSNTOP(JL)>=PTSN(JL)) THEN
        ZACTDEPTH(JL)=ZSNDEPTH(JL, KLEVSNA(JL))
      ELSE
        ZACTDEPTH(JL)=ZDSNTOT(JL)
      ENDIF

        PTSNMIDDLE(JL)= MIN(RTT,0.5_JPRB*(PTSN(JL)+PTSOIL(JL)))
        PTSNBOTTOM(JL)= MIN(PTSOIL(JL),RTT)
!----- Density
      PRSNTOP(JL)=300._JPRB
    ENDIF
  ENDDO 

END ASSOCIATE

!    -----------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('SURFWS_INIT_SL_MOD:SURFWS_INIT_SL',1,ZHOOK_HANDLE)

END SUBROUTINE SURFWS_INIT_SL
END MODULE SURFWS_INIT_SL_MOD


